Applications of SOS Programming in Flowpipe Construction

Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated XFZ18under). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.

In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using XFZ18under. This is a Julia implementation that relies on the JuMP ecosystem (JuMP, PolyJuMP, SumOfSquares, MathProgInterface, MathOptInterfaceMosek) and the JuliaAlgebra ecosystem (MultivariatePolynomials, DynamicPolynomials). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.


References:

Quick links to documentation:

Van-der-Pol system

Consider the following van-der-Pol system:

$$ \dot{x}_1 = x_2 \\ \dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2 $$

In [1]:
using MultivariatePolynomials,
      JuMP,
      PolyJuMP,
      SumOfSquares,
      DynamicPolynomials,
      MathOptInterfaceMosek,
      MathematicalSystems

const  = differentiate
Out[1]:
differentiate (generic function with 18 methods)
In [79]:
# symbolic variables
@polyvar x₁ x₂ t

# time duration (scaled, see dynamics below)
T = 1.0 

# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂] 

# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25

# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 9 - x₁^2 - x₂^2

# degree of the relaxation
k = 12

# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)

# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(MosekOptimizer))

# add unknown Φ to the model
@variable(model, Φ, Poly(XT))

# jacobian
∂t = α -> (α, t)
∂xf = α -> (α, x₁) * f[1] + (α, x₂) * f[2] 
 = ∂t(Φ) + ∂xf(Φ)

# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)

# scalar variable
@variable(model, ϵ)

dom1 = @set t*(T-t) >= 0 && g >= 0
dom2 = @set g >= 0
@constraint(model, ϵ >= 0.)
@constraint(model,   SOSCone(), domain = dom1)
@constraint(model, ϵ -   SOSCone(), domain = dom1)
@constraint(model, Φ₀ - V₀  SOSCone(), domain = dom2)
@constraint(model, ϵ + V₀ - Φ₀  SOSCone(), domain = dom2)

@objective(model, Min, ϵ)
Out[79]:
$$ ϵ $$
In [80]:
optimize!(model)
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31448) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31449) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31450) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31451) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31452) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31453) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31454) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31455) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31456) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(31457) of matrix 'A'.
Warning number 705 is disabled.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 31617           
  Cones                  : 0               
  Scalar variables       : 30530           
  Matrix variables       : 10              
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 396
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.01            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.12    
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 31617           
  Cones                  : 0               
  Scalar variables       : 30530           
  Matrix variables       : 10              
  Integer variables      : 0               

Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 31217
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 30132             conic                  : 30131           
Optimizer  - Semi-definite variables: 10                scalarized             : 30074           
Factor     - setup time             : 69.29             dense det. time        : 8.29            
Factor     - ML order time          : 6.44              GP order time          : 37.62           
Factor     - nonzeros before factor : 7.85e+07          after factor           : 9.63e+07        
Factor     - dense dim.             : 2                 flops                  : 4.18e+11        
Factor     - GP saved nzs           : 1.37e+08          GP saved flops         : 1.73e+12        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  4.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  69.44 
1   5.9e-01  2.4e+00  8.8e-01  1.47e+00   6.535729443e-01   1.366129999e-02   5.9e-01  78.24 
2   1.3e-01  5.4e-01  9.0e-01  1.85e+00   5.505220208e-01   2.681089005e-02   1.3e-01  87.17 
3   2.5e-02  1.0e-01  8.4e-01  2.05e+00   1.564279934e-01   8.659716489e-02   2.5e-02  96.48 
4   2.5e-03  1.0e-02  4.4e-01  1.50e+00   2.182888976e-02   1.680843274e-02   2.5e-03  109.12
5   3.3e-04  1.3e-03  1.5e-01  1.25e+00   6.396295208e-03   5.805212257e-03   3.3e-04  119.49
6   4.4e-05  1.8e-04  4.2e-02  9.51e-01   2.074509936e-03   1.992919035e-03   4.4e-05  129.88
7   3.4e-05  1.4e-04  3.5e-02  6.72e-01   1.805736560e-03   1.736794232e-03   3.4e-05  140.76
8   5.8e-06  2.3e-05  9.5e-03  6.97e-01   1.053338983e-03   1.039659671e-03   5.8e-06  151.08
9   3.3e-06  1.3e-05  4.5e-03  8.92e-03   1.586165278e-03   1.575837809e-03   3.2e-06  160.27
10  8.4e-07  2.7e-06  7.9e-04  -1.07e-01  3.188592768e-03   3.190188702e-03   6.7e-07  170.20
11  6.2e-07  1.6e-06  4.2e-04  -3.09e-01  4.269761781e-03   4.275329009e-03   3.9e-07  179.10
12  1.7e-07  4.4e-07  9.2e-05  -3.18e-01  7.679991847e-03   7.696992349e-03   1.1e-07  200.41
13  1.0e-07  1.6e-07  3.5e-05  3.10e-02   1.114944880e-02   1.116667699e-02   3.9e-08  219.06
14  1.0e-07  1.6e-07  3.5e-05  9.74e-02   1.114944880e-02   1.116667699e-02   3.9e-08  230.22
Optimizer terminated. Time: 231.69  

In [81]:
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
Relaxation order : k = 12
JuMP.termination_status(model) = SlowProgress
JuMP.primal_status(model) = FeasiblePoint
JuMP.dual_status(model) = FeasiblePoint
JuMP.objective_bound(model) = 0.0
JuMP.objective_value(model) = 0.011149448804550344
In [45]:
# Recovering the solution:
ϵopt = JuMP.objective_value(model)

# Punder <= 0
Punder = subs(JuMP.value(model[:Φ]), t => T)
Out[45]:
-1.2042160063789687e-5x₁^{10} - 4.2399855141005857e-5x₁^{9}x₂ + 5.397739071992859e-5x₁^{8}x₂^{2} - 0.0003281664096657183x₁^{7}x₂^{3} + 0.00027514830691699845x₁^{6}x₂^{4} - 0.0004299553713378092x₁^{5}x₂^{5} + 0.00022286405919632663x₁^{4}x₂^{6} - 0.00018873068043793128x₁^{3}x₂^{7} + 6.950578607376682e-5x₁^{2}x₂^{8} - 2.8081627666802042e-5x₁x₂^{9} + 2.9498267276158296e-6x₂^{10} + 4.471679594848281e-20x₁^{9} - 6.38719542038995e-20x₁^{8}x₂ + 1.5280606303157782e-19x₁^{7}x₂^{2} + 4.233615217434361e-20x₁^{6}x₂^{3} - 3.604204090946449e-19x₁^{5}x₂^{4} + 3.4899346915783347e-19x₁^{4}x₂^{5} - 2.9856665228720004e-19x₁^{3}x₂^{6} + 1.8772310102595033e-19x₁^{2}x₂^{7} - 9.444126451164612e-20x₁x₂^{8} + 1.8406027322704353e-20x₂^{9} + 0.001895650996724096x₁^{8} - 0.0013293724244520244x₁^{7}x₂ - 0.001948586395811546x₁^{6}x₂^{2} + 0.0075236402265226494x₁^{5}x₂^{3} - 0.00507770572933377x₁^{4}x₂^{4} + 0.005395129670341276x₁^{3}x₂^{5} - 0.002666130050900258x₁^{2}x₂^{6} + 0.0011262027500646237x₁x₂^{7} - 0.00012479887873173822x₂^{8} + 1.2351542659507538e-18x₁^{7} - 2.713031271959235e-18x₁^{6}x₂ + 2.390246389072238e-18x₁^{5}x₂^{2} - 2.228238381291506e-18x₁^{4}x₂^{3} + 4.437390326050367e-18x₁^{3}x₂^{4} - 4.4432504374509434e-18x₁^{2}x₂^{5} + 2.4312072694839196e-18x₁x₂^{6} - 4.964935509264464e-19x₂^{7} - 0.03545590192538548x₁^{6} + 0.06545571953996324x₁^{5}x₂ - 0.03567152236694011x₁^{4}x₂^{2} - 0.03697234276328999x₁^{3}x₂^{3} + 0.034079735366444305x₁^{2}x₂^{4} - 0.015031050502646914x₁x₂^{5} + 0.00202443363386639x₂^{6} - 2.153249801569398e-17x₁^{5} + 4.500204544676722e-17x₁^{4}x₂ - 5.674641884575686e-17x₁^{3}x₂^{2} + 4.28076703948008e-17x₁^{2}x₂^{3} - 2.0652882452373237e-17x₁x₂^{4} + 3.979192608630357e-18x₂^{5} + 0.19253582548215603x₁^{4} - 0.5529398295468485x₁^{3}x₂ + 0.7216844609906954x₁^{2}x₂^{2} - 0.3867031818373412x₁x₂^{3} + 0.08587589025436114x₂^{4} + 7.177602168893349e-17x₁^{3} - 8.305255834222396e-17x₁^{2}x₂ + 4.576059730878753e-17x₁x₂^{2} - 7.20361079584399e-18x₂^{3} + 0.1703861869394433x₁^{2} - 0.518073359580279x₁x₂ + 0.29621187452208164x₂^{2} - 5.260357888653161e-17x₁ + 1.0615570214913224e-17x₂ + 0.1330556718203539
In [46]:
# Pover <= 0
Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)
Out[46]:
-1.2042160063789687e-5x₁^{10} - 4.2399855141005857e-5x₁^{9}x₂ + 5.397739071992859e-5x₁^{8}x₂^{2} - 0.0003281664096657183x₁^{7}x₂^{3} + 0.00027514830691699845x₁^{6}x₂^{4} - 0.0004299553713378092x₁^{5}x₂^{5} + 0.00022286405919632663x₁^{4}x₂^{6} - 0.00018873068043793128x₁^{3}x₂^{7} + 6.950578607376682e-5x₁^{2}x₂^{8} - 2.8081627666802042e-5x₁x₂^{9} + 2.9498267276158296e-6x₂^{10} + 4.471679594848281e-20x₁^{9} - 6.38719542038995e-20x₁^{8}x₂ + 1.5280606303157782e-19x₁^{7}x₂^{2} + 4.233615217434361e-20x₁^{6}x₂^{3} - 3.604204090946449e-19x₁^{5}x₂^{4} + 3.4899346915783347e-19x₁^{4}x₂^{5} - 2.9856665228720004e-19x₁^{3}x₂^{6} + 1.8772310102595033e-19x₁^{2}x₂^{7} - 9.444126451164612e-20x₁x₂^{8} + 1.8406027322704353e-20x₂^{9} + 0.001895650996724096x₁^{8} - 0.0013293724244520244x₁^{7}x₂ - 0.001948586395811546x₁^{6}x₂^{2} + 0.0075236402265226494x₁^{5}x₂^{3} - 0.00507770572933377x₁^{4}x₂^{4} + 0.005395129670341276x₁^{3}x₂^{5} - 0.002666130050900258x₁^{2}x₂^{6} + 0.0011262027500646237x₁x₂^{7} - 0.00012479887873173822x₂^{8} + 1.2351542659507538e-18x₁^{7} - 2.713031271959235e-18x₁^{6}x₂ + 2.390246389072238e-18x₁^{5}x₂^{2} - 2.228238381291506e-18x₁^{4}x₂^{3} + 4.437390326050367e-18x₁^{3}x₂^{4} - 4.4432504374509434e-18x₁^{2}x₂^{5} + 2.4312072694839196e-18x₁x₂^{6} - 4.964935509264464e-19x₂^{7} - 0.03545590192538548x₁^{6} + 0.06545571953996324x₁^{5}x₂ - 0.03567152236694011x₁^{4}x₂^{2} - 0.03697234276328999x₁^{3}x₂^{3} + 0.034079735366444305x₁^{2}x₂^{4} - 0.015031050502646914x₁x₂^{5} + 0.00202443363386639x₂^{6} - 2.153249801569398e-17x₁^{5} + 4.500204544676722e-17x₁^{4}x₂ - 5.674641884575686e-17x₁^{3}x₂^{2} + 4.28076703948008e-17x₁^{2}x₂^{3} - 2.0652882452373237e-17x₁x₂^{4} + 3.979192608630357e-18x₂^{5} + 0.19253582548215603x₁^{4} - 0.5529398295468485x₁^{3}x₂ + 0.7216844609906954x₁^{2}x₂^{2} - 0.3867031818373412x₁x₂^{3} + 0.08587589025436114x₂^{4} + 7.177602168893349e-17x₁^{3} - 8.305255834222396e-17x₁^{2}x₂ + 4.576059730878753e-17x₁x₂^{2} - 7.20361079584399e-18x₂^{3} + 0.1703861869394433x₁^{2} - 0.518073359580279x₁x₂ + 0.29621187452208164x₂^{2} - 5.260357888653161e-17x₁ + 1.0615570214913224e-17x₂ - 0.2855694195323183
In [47]:
using ImplicitEquations, Plots
gr() # or pyplot()

G = plot()

_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
plot!(G, _Punder  0., xlims=(-10, 10), ylims=(-10, 10), color="red")

_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
plot!(G, _Pover  0., xlims=(-10, 10), ylims=(-10, 10), color="blue")

G
Out[47]:
-10 -5 0 5 10 -10 -5 0 5 10
In [25]:
G = plot()

_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
Gu = plot(_Punder  0., xlims=(-8, 8), ylims=(-8, 8))

_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
Go = plot(_Pover  0, xlims=(-8, 8), ylims=(-8, 8))

plot(Gu, Go)
Out[25]:
-6 -3 0 3 6 -6 -3 0 3 6 -6 -3 0 3 6 -6 -3 0 3 6

Comparison with a MATLAB/YALMIP implementation

Package k Constraints Scalar variables Matrix variables Time(s)
JuMP 2 245 175 8 < 1
YALMIP 2 152 63 10 < 1
JuMP 3 893 715 10 ~ 1
YALMIP 3 254 121 10 ~ 1
JuMP 4 893 730 10 ~1
YALMIP 4 394 206 10 1.18
JuMP 5 2639 2309 10 1.17
YALMIP 5 578 323 10 0.11
JuMP 6 2639 2337 10 0.90
YALMIP 6 812 477 10 1.10
JuMP 7 6725 6183 10 5.62
YALMIP 7 1102 673 10 1.52
JuMP 8 6725 6228 10 6.06
YALMIP 8 1454 916 10 1.10
JuMP 9 15269 14447 10 41.2
YALMIP 9 1874 1211 10 1.58
JuMP 10 15269 14513 10 38.1
YALMIP 10 2368 1563 10 2.73
JuMP 11 31617 30439 10 244
YALMIP 11 2942 1977 10 2.30
JuMP 12 31617 30530 10 231
YALMIP 12 3602 2458 10 6.57